Biaxial spin-nematic phase of two dimensional disordered rotor models and spin-one bosons in 

optical lattices 
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We show that the ground state of disordered rotor models with quadrupolar interactions can exhibit biaxial 
nematic ordering in the disorder-averaged sense. We present a mean-field analysis of the model and demonstrate 
that the biaxial phase is stable against small quantum fluctuations. We point out the possibility of experimental 
realization of such rotor models using ultracold spin-one Bose atoms in a spin-dependent and disordered optical 
lattice in the limit of a large number of atoms per site and also suggest an imaging experiment to detect the 
biaxial nematicity in such systems. 
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I. INTRODUCTION 

Nematic phases of matter have been widely studied in many 
different areas of condensed matter physics^. For example, 
consider a system of molecules having some vector quantity 
(say a polarization vector p) associated with each of them. 
Since molecules are extended objects, a local coordinate sys- 
tem can be attached to each molecule. Hence, for a molecule 
at a point i, one can define a local nematic order parame- 
ter: = ViaV%0 - HaVlJ^-- The nematic phase is 
then defined as the state of the system for which £\ p; = 

whereas J2i Qi becomes a traceless 3x3 matrix with eigen- 
values Q = [—(Qi + Q2), Qi, Q2Y Such phases are of two 
types: uniaxial nematic where Q\ = Q2 and biaxial ne- 
matic for which Q\ ^ Q2 1 ■ In both these phases, the sys- 
tem breaks rotational symmetries while preserving transla- 
tional invariance. Uniaxial nematic phases occurs more fre- 
quently whereas the first experimental observation of biax- 
ial nematic phase in a classical system of molecules has only 
been recently reported 2 . Although much more difficult to re- 
alize, biaxial nematic phases are of great theoretical interest 
since they are known to support non-Abelian defects 3 . Uni- 
axial nematic phases have also been studied in the context of 
several quantum condensed matter systems such as strongly 
correlated electron system a 4 ' 5 ' 6 ' 7 ' 8 ' 9 ' 10 and ultracold spin-one 
atoms in optical lattic e) 11,12 ' 13 . For example, in the case of 
spin-one boson system, in some parameter regime, the ground 
state breaks spin-rotational invariance along one of the axes 
while maintaining translational invariance with no spin-order 
(S) = and with two equal non-zero eigenvalues of the ma- 
trix (S a SP). However, there has been no proposals of realiza- 
tion of biaxial nematic ground states in the context of quantum 
condensed matter systems using standard local Hamiltonians. 

In this work, we point out that disordered O (2) rotor models 
with quadrupolar interaction in two dimensions (2D) can ex- 
hibit biaxial nematic phase in the disordered averaged sense. 
Such rotor models are described by the Hamiltonian 



A a d- 



(1) 



<ij> \<^—x,y,z 



where d = (d x ,d y ,d z ) = (sin 6 cos <f), sin 6 sin <f>, cos ff) is a 
unit vector expressed in terms of the rotor angles and <f>, 
Li denotes angular momentum of the rotor, Yluj) denotes 
sum over nearest neighbor sites, and Ay- are dimensionless 
random couplings between the unit vectors d at sites i and 
j. In the following we shall consider A x = A y ^ A z , and 
choose A x ^ to have Gaussian distributions with same width 
SA and means A x ^ . The main result of this work is that in the 
limit \A X — A z \ /SA <C 1, the model given by Eq.[TJ exhibits 
a biaxial nematic ground state in the disordered average sense 
for U -C 5- We also show, using a mean-field analysis, that 
the biaxial phase is stable against quantum fluctuations until 
v = g/U is reduced to a critical value v c . Note that the above- 
mentioned regime can also be reached for arbitrarily weak dis- 
order: SA < A X ,A Z provided that \A X - A z \ < SA. This 
particular feature of the rotor model, as we shall see, makes it 
easily realizable using ultracold atoms in optical lattices. 

The choice of the rotor Hamiltonian (Eq.[TJ) is not arbitrary, 
but is motivated by their connection with ultracold spin-one 
bosons in an optical lattice. The low energy effective Hamil- 
tonian of such bosons systems has been derived in the ab- 
sence of disorder in Refs. [T11ll2lll3L It has been shown that 
the effective theory of the Mott phases of such a system can 
be described, in the limit of large number of bosons per site, 
by rotor models similar to Eq. [TJ supplemented by an addi- 
tional constraint on the number and the total spin S of the 
bosons per site: N + S = eve n 11 ' 12 ' 13 . We generalize this 
analysis to include disorder and show that in the presence of 
spin-dependent disordered optical lattices, which can be easily 
created by using a speckled laser fiel d 14 ' 15 , the Mott phases of 
ultracold spin-one bosons is indeed described by Eq.[TJsupple- 
mented by the constraint condition mentioned above. Further, 
as will be shown later, for U <C g, the constraint condition 
becomes irrelevant and in this limit we expect our analysis of 
the rotor model to apply for the bosons as well. This corre- 
spondence therefore allows us to suggest experiments on ul- 
tracold atoms which can in principle probe such a biaxial ne- 
matic phase. We suggest a straightforward imaging of ultra- 
cold atoms using a polarized laser beam which can distinguish 
the biaxial phase from other spin or uniaxial nematic phases. 
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The organization of the paper is as follows. In Sec. [II] we 
analyze the rotor model (Eq. [TJ and demonstrate that it ex- 
hibits a biaxial nematic ground state. This is followed by Sec. 
IIII Al where we derive the effective low-energy Hamiltonian 
of spin-one ultracold Bosons in a spin-dependent disordered 
optical lattice and discuss its relation to the rotor model an- 
alyzed in Sec. [ID In Sec lIIIBI we suggest an experiment to 
detect biaxial nematic order. Finally, we conclude in SeclIVI 

II. ANALYSIS OF THE ROTOR MODEL 

In this section, we first present an analysis of the rotor 
model (Eq. [TJ and demonstrate the presence of biaxial ne- 
maticity at U = 0. This is followed by the derivation and 
analysis of a mean field theory for finite U which illustrates 
the stability of the biaxial phase against small quantum fluc- 
tuations. 



A. Limit of U = 

For U = 0, in the absence of quantum fluctuations, the 
Hamiltonian (Eq,[T]i is diagonal in the \Q, 4>) basis. This allows 
us to write 

where Ay is given by 

Aij = A^j sin 0i sin 9j cos(0i — <f>j) + A\j cos Qi cos 9j (3) 

with < 8i < 7r and < <j>i < 2tt. Let us now consider a 
ground state configuration of the Hamiltonian H' mtor given by 

* G = \e\, <t>\;6l, 0°; 0°; ....6%, </>%• ....) (4) 

where (#°, 0°) denotes the value of the angles Q and </> at the 
ith site in the ground state. From Eqs.[2]and[3] one can imme- 
diately see that H' rotor has an infinite number of degenerate 
ground states since any local transformation 8® — > it — 6®, 
(j)® — > <f>1 + 7T leaves H[ otm invariant for any random set of 
Af; 2 . It is also easy to see that such a local transformation 
changes — > d; and leaves the local nematic order param- 
eter Qf b — dfcfy — Sab/3 invariant. Hence for the purpose 
of computing the nematic order parameter, we can choose any 
one of these ground state configurations. 

The rotor Hamiltonian H' rot is also invariant under two 
global transformations. The first of them T\ : <f> — * <f> + 77 is 
a gauge freedom which allows us to choose the orientation of 
the global x axis. The second transformation T 2 : Q — > it — Q 
leaves the diagonal components of Q l ab invariant while chang- 
ing the off-diagonal components. The ground states ^g (Eq. 
D and T 2 ^g = *' G ^ degenerate. 

Using the above-mentioned local and global transformation 
properties, we can considerably simplify the numerical proce- 
dure for obtaining ground states of H' TotoT for a given disorder 
realization. First, from Eq. [3] we find that any ground-state 



configuration $g will always have <f)i — <f)j — 0, ir. Using the 
transformation T\ , we can therefore choose each individual 
cf>i to be either or ir. This is a gauge choice and, hence, does 
not affect any physical quantity. Next instead of using ^>g to 
compute Q ab , we find the state for which <fii — at every 
site. This can be done since \& G is always a part of the degen- 
erate states which can be reached from ^g by successive ap- 
plication of the above-mentioned local transformation. Using 
these two facts, we see that for any set of random coefficients 
A^ 2 , one can compute Q ab = J2i A la d ih - 6 ab /3 by finding 

the ground state configuration of H TotoI = -gJ2<ij>(^'ij) 2 
with 

A'^ = sin 8, sin 9j + A^ cos 6»; cos 6j . (5) 

The ground states of H mtoI is two-fold degenerate and these 

two ground states \& G ^ and = T2^g mQ related by 
the global transformation T 2 . The two ground states have 
opposite signs for off-diagonal elements of Q a b in the cho- 
sen gauge and therefore after averaging over these ground 
states, the off-diagonal elements of Q ao vanishes. Thus 
we are finally left with the diagonal components Q aa = 
(-(Qi + Q 2 ),Qi,Q 2 ) - [C- 1/3, -1/3, -(C- 2/3)), 
where < C < 1, for a given set of random coefficients 
A^' 2 . Finally, one can numerically average over different re- 
alizations Q a = (Qaa)disorder t0 obtain the disorder-averaged 
nematic order parameter. The biaxial nematic phase is real- 
ized when C ^ 0, \ or 1. 

The ground state of H" 10 t r can be obtained numerically 
using standard minimization procedure for each disorder real- 

x( z) 

ization. For the sake of brevity, we choose A^ to be a set 
of random numbers having a Gaussian distribution with same 
standard deviation SA and different means A x ^ and average 
over 500 disorder realizations. We choose A^ > A z here 
without loss of generality. We study finite-size systems with 
sizes ranging from N$ — 8 x 8 to N$ = 20x20 with periodic 
boundary conditions and use standard 1 / N$ extrapolation to 
obtain the infinite system size limit 16 . A few comments about 
the numerical procedure are in order here. 

First, for (A x - A Z )/SA > 1, it is clear that the sys- 
tem energy reaches its minimum at sin Qi sin 6j = 1 and 
cos Qi cos Qj — for all sites and disorder realizations. Conse- 
quently, after disorder-averaging, (7=1, and the ground state 
is an uniaxial nematic phase. 

Second, in the limit (A x — A Z )/6A 1, for a given disor- 
der realization, we find that the configuration with minimum 
energy is part of a group of three possible solutions. The 
first corresponds to sin Qi sin Qj = 1 and cos Qi cos Qj = 
for all sites. For a single disorder realization with this solu- 
tion, the order parameter then has C = 1. The second solu- 
tion possibly assumed by the system is sin Qi sin Qj = and 
cos Qi cos Qj ~ ±1 for all sites. The value of the order param- 
eter corresponding to such a disorder realization is then given 
by C = 0. Finally, the third possible solution corresponds to a 
configuration for which sin Qi sin Qj and cos Qi cos Qj are nei- 
ther nor ±1, and are varying in values from site to site. For 
such a disorder realization, C 7^ 0, ^ or 1 and one has a biaxial 
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(A x — A Z )/8A, the biaxial phase will manifest itself over a 
large window of energy scale at finite temperature, even if the 
zero temperature ground state turns out to be uniaxial nematic. 



B. Quantum fluctuations for U 7^ 

To study the effect of quantum fluctuations, it is more con- 
venient to switch to a path integral representation of the rotor 
Hamiltonian (Eq. [T). To this end, we first write the partition 
function of our Hamiltonian as 
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FIG. 1: — (Qi + Q2) = C — 1/3 as a function of the inverse lattice 
size for several disorder strengths. For moderate disorder (A x — 
A Z )/5A > 0.1, an extrapolation to large system size shows that the 
ground-state will be uniaxially ordered, C — ► 1. For strong disorder 
(A — A Z )/SA < 0.07, the ground-state remains biaxially ordered, 
C / 1 after disorder averaging within 1/Ns extrapolation. [Note: 
Qi = <Qi}, Q2 = (Q2), C = (C), A x = (A*}, A z = <A Z }] 



nematic. However, having C 7^ 0, i or 1 for a single disor- 
der realization does not guarantee that the system is biaxially 
ordered. The system can in fact be separated into domains, 
each domain representing one of the three possible solutions. 
For example, the system could be separated into domains rep- 
resenting only the first and second type of solutions and still 
have C ^ 0, I or 1. Nonetheless, this situation is unlikely; 
in most cases where C 7^ 0, \ or 1 we find that the system 
adopts the third solution. Nevertheless, as mentioned previ- 
ously, to find the physically meaningful ground state, we need 
to average over many disorder realizations and by doing so re- 
store the translational invariance of the system. Consequently, 
C 7^ 0, \ or 1 can be used only after disorder-averaging to 
determine whether the system is biaxially ordered for certain 
disorder strength. 

The result of the numerical calculation, which corrobo- 
rates the above discussion, is shown in Fig. Q] After fitting 
second order polynomial functions to our data and extrap- 
olating the resulting functions to 1/Ns — > 0, we find that 
for (A x — A Z )/SA < 0.07, the ground state is biaxial ne- 
matic with C / 1. For intermediate or larger values of 
(A x — A Z )/5A > 0.1, C = 1 and we have an uniaxial ne- 
matic ground state. So this demonstrates that for U = 0, the 
rotor model exhibits a biaxial nematic ground state in the dis- 
ordered average sense in the limit (A x — A Z )/SA <C 1. In the 
next subsection, we study the effect of quantum fluctuations 
on this state. 

A comment about the validity of the 1/Ns extrapolation is 
in order here. This procedure definitely leaves open a possibil- 
ity that a large enough system size has not been reached in our 
numerics and the biaxial nematic state seen here is an artifact 
of finite system size. The resolution of this question, within 
our numerical procedure, is not straightforward. However, we 
can be certain that even if this is the case, with small enough 



Z = Tr cxp(-/3H rot or) = Tr exp [—(3(T + V)} (6) 
T = Uj2 L i (7) 



abed (ij) 



(8) 



where we have introduced the notation = Af-Af-S a hS cc i. 
To obtain the partition function we calculate the trace in (|6]l 
by writing it as a path integral over M time slices, and then 
inserting a complete set of coherent states \6,<j)) over each 
slice, so that 



Z = lim Tr[exp{-(5r(T + F)}] 

M— »oo 

M-l 



M 



jv9V4> Yl (O(T a+1 );0(T a+1 )\cxp[-8TT] 
x exp [StV] |0(t o ); 0(r a )) 



(9) 



where St = /3/M. Inserting a set of complete states \l,m) 
(defined by (8, <j>\l, m) = Yi m (8, </>), where Yi m are the spher- 
ical harmonics) on each slice, we get 



M-l 



V8V(j) Y[ exp [StV (8, </>)] T a 



(10) 



T a = J2 Jl Hexp[-USTl a (i)(l a (i) + l)} 



l a =0 7 



--la, i 



X Y C (i),m a (i) tfa+l {i),4> a +l{i)) 

xYi a (i) tm , a (i)(0 a (i),(l) a (i)) (11) 

To rewrite (fTTT i in a suitable form, we use the identities 17 

00 / 



s(A0) 



= VtE E I lH {h)YC m {8^)Y lm {0'^') 



1=0 m=— I 



lim I l+ i{h) = exp 



(/ + l/2) 2 -l/4 
2h 



0(l/h 2 ) 



(12) 



where I l+ i is the modified Bessel function and cosA6> = 
cos 8 cos 8' + sin sin 0' cos(0 — 4>'). Then, after a few 
straightforward manipulations, we get in terms of d = 



4 



(sin 9 cos <fi, sin 9 sin 0, cos 9) fields, 

Z = J Dd^S (^2\d ia \ 2 - lj exp(-S) 



(13) 



rfr 



5^ (a T d w ) 2 - v E d ia d tc r% c % b d Jd 

ia <ij> abed 

(14) 

where we have rescaled St so that 4U St = 1 and v = g/ 4U. 

Next, we introduce a Lagrange multiplier field A,; (r) to im- 
plement the constraint and decouple the quartic term in S us- 
ing an auxiliary field N l ab . After some algebra, we get 



Z = J T>dT>NV\exp(-Si) 



(15) 



Si = 



/ dT E £(0rd») 2 + iAi( ^|rf io | 2 -lj 

L i L a \ a ) 



a b 



<ij> abed 



(16) 

where the auxiliary fields AT* 6 are not the order parameter 
fields, but their conjugate. The order parameter fields P % ab can 
now be introduced by a second Hubbard-Stratonovitch trans- 
formation which yields 



Z = J VdVNVXVPcxp{~S cfi ) 



(17) 



S, 



off 



J dT E { E ( 9 - d -) 2 + a > (e i d -i 2 - : ) 



<ij> abed 



(18) 



An integration over the auxiliary fields N* b shows that P l ab = 
(diadib) and hence the nematic order parameter Q ab — 
Si ( P ab ~ $ab J2 C P cc/3>) can be directly obtained in terms 
of the P* b fields. 

Consequently, we can now seek the saddle point solution 
to the above action. At the saddle point, the constraint fields 
are time-independent and the d; fields can be integrated out. 
Notice that in contrast to the clean system, the constraint fields 
\i are space-dependent. The mean-field action becomes 



Smf = E Tr [ lnG W] 



J dT j2(^i-iJ2 N ab p i 

L i ab 



(19) 



<ij> abed 

[{-d 2 T + S ab ~ iN l ab ] (20) 




100 



FIG. 2: -(Qi+Q 2 )asafunctionoiV = g/AU for (A x -A Z )/5A = 
0.03. This result, obtained for N s = 10 x 10 and N s = 12 x 12 
lattices, shows that for a wide range of v > v c the ground-state 
remains biaxial nematic. [Note: Q\ = (Qi), Qi = (Q2)] 



The saddle point equations can now be obtained from Eq.[T9l 

1 = /^E[ G f/M^)] (2D 



iN, 



ab 



r ab — 



V E E ^acbd P cd 
cd j 

du> 



2n 



(22) 
(23) 



and are solved numerically to obtain the mean field order pa- 
rameter Q ab = J^i ( P ab ~ S a,b J2 C P ld?>) for each disorder 
realization. 

For a given disorder realization, we solve the set of mean- 
field equations Eqs.|2Tl|22j and [23] in order to find, for each 
lattice site i, the fields P* b (with a.b = {x, y, z}). We then 
average over 100 disorder realizations and for finite sizes 
N s = 10 x 10 and 12 x 12. We find that the real val- 
ued solutions of the order parameter are diagonal and obtain 

Qa = (Qa)di SO rder = +Qi),Qi,Qi)- The result is 

shown in Fig. for (A x - A Z )/SA = 0.03. We find, from the 
plot of — (Qi + Q2) as a function of v, that the biaxial phase 
persists for a wide range of v > v c = 1. This shows that the 
effect of quantum fluctuations, at least at a saddle point level, 
does not destabilize the biaxial phase as long as v > v c . This 
qualitative feature seems to be independent of system size, 
as can be seen from Fig. [2] and we expect it to hold in the 
Ns — > 00 limit. 



III. ULTRACOLD SPIN-ONE ATOMS 

In this section we first show that the low energy effective 
Hamiltonian of spin-one ultracold atoms in the Mott insulat- 
ing phase can be mapped onto the rotor Hamiltonian (Eq.Q]) 
which we analyzed earlier. Then we propose an imaging ex- 
periment on ultracold atoms which can detect the biaxial ne- 
matic phase. 
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A. Effective rotor model 

We consider a system of bosonic 5=1 atoms in a disor- 
dered optical lattice with spin-dependent confining potentials 
and antiferromagnetic interaction between atoms. The many- 
body Hamiltonian for this system is given in second-quantized 
notation by 



H 



dr &(r) (-^-W 2 + V a {v)\ &(r) 

+ \ J drdr'$l(r)$l,(r')W(r-r')Mr')i> b ,(r), 

(24) 

where ip\ (r) is the boson field operator that creates a particle 
with spin projection a = { — 1,0, 1} at position r, V a (r) is 
the spin-dependent and spatially disordered external potential, 
and W(r - r') = <5(r - r')(U + U 2 Sx • S 2 ) is the two- 
body interatomic potential 1 ^. Here £/ 2 = ^^~( a 2 — a-o) and 

Uq = ^^-(2a2 + a ) are the on-site interactions, a 2 (o) are 
the s -wave scattering lengths in the S = 2(0) spin channels, 
and m is the mass of the atoms. For example, for 23 Na, the 
scattering lengths are a 2 = (52 ± 5)as and clq = (46 ± 5)as 
where as is the Bohr radius 12 so that U2/U0 ~ 0.04. In what 
follows, we shall be interested in the Mott states of the spin- 
one bosons which occur in the limit of deep lattice potential 
with Vx = V-x £ V . 

A spin-dependent disordered potential V Q (r) can be gen- 
erated by superposing a speckled laser field on a sinusoidal 
spin-dependent lattice potential Vba(r). The spin dependence 
of the sinusoidal potential can be achieved by tuning the laser 
frequency close to the hyperfine splitting (but far away from 
the fine structure splitting) of the atoms 19 i 20 i 2 1 ,22.23.24 ^ j n a p_ 
pendix[A] we propose a method to obtain a spin-dependent 
lattice potential with V\ = V-\ 7^ Vo- Generation of the dis- 
ordered potential 5 V a (r) can be achieved by reflecting a laser 
at the same frequency but with a much lower intensity off a 
speckled mirror. The net lattice potential seen by the atoms at 
spin state a is then V a (r) — Vo a (r) + SV a (r). In what fol- 
lows, we shall consider the situation where the speckled field 
laser is weak compared to the one generating the sinusoidal 
potential such that 5V a <C Vo a - 

For free ultracold atoms in an optical lattice, the energy 
eigenstates are Bloch wave functions and a superposition of 
these Bloch states yields a set of Wannier functions which are 
well localized on the individual lattice sites for deep lattices 25 . 
The energies involved in the system dynamics being small 
compared to excitation energies to the second band, we ex- 
pand the boson field operators in the Wannier basis and keep 
only the lowest states, V'a(r) = Yli bi a w(r — r^). Conse- 
quently, the many-body Hamiltonian d24T i reduces toii 



where bi a is the spin-one boson operator at site i with spin pro- 
jection a = { — 1, 0, 1}, hi = J2 a "iaf>ia is tne boson density 
at site i, S, = ^ b b\ a S a bbib is the spin operator (S being the 
spin rotation matrices for spin-one bosons), fj, is the chemical 
potential, and are given by 



drw*(r-r t ) ( + K(r) ) ^(r-r,). (26) 



so that tx = t-x ^ h for Vx = V-x ^ V . 

Note that for weak speckled fields, we always have 
o't/ia "C 1, where at and t a are the standard deviation and 
average of the hopping coefficient t a . However, in this setup, 
one can always tune the sinusoidal potential Vo a so that we 
are in a regime where \ix — <o|/ cr t *C 1- As we shall see, this 
is precisely the regime where one expects to see the biaxial 
phase. Also we note that due to the on-site disordered poten- 
tial, the chemical potential ^ should also be site-dependent. 
However, //j has only a power law dependence on disorder by 
comparison to an exponential dependence for the hopping co- 
efficients. Consequently, the standard deviation of fa is small 
and does not influence the nature of the Mott states. So we re- 
placed fa by its average value /i in Eq.|25] Notice that this ap- 
proximation is valid only when the potential due to the speck- 
led field is weak compared to the one due to the sinusoidal 
field. 

Next, following Ref. [TH we switch to a representation 
where the boson operators transform as vectors under spin ro- 
tation 

1 —i 
hz = b i0 , b ix = — (6, _i - bix), b iy = — (6i_i + hx). 

(27) 

Using these operators, the kinetic term in (|25] l is rewritten as 

(28) 



(ij) a£{x,y,z} 



with t'J = = 2t[ j and t\i = itf. 

We now consider this spin-one Bose-Hubbard model in the 
limit where N ^> 1 and we are in the Mott phase of the bosons 
with a large odd number (N) bosons per site. A straightfor- 
ward generalization of the analysis of Ref. [TH shows that the 
low energy effective Hamiltonian in this limit can be mapped 
on to a rotor model. Using the decomposition b La = di a di, 
where the boson operator dj changes the number of particles 
Ni but not the orientation of the boson spin given by d;, we 
obtain in second order perturbation theor y 11 ' 12 



K 



off 



— D ( T, ^h- 

<ij> \a=x,y,z x / 



U 



(29) 



(25) 



(ij) a 



Eq.|29]has to be supplemented with the constraint Ni + Si = 
even where Ni is the number of bosons at site i, and Si — 
ieijkdj gf- is the total spin which can be identified as the rotor 
angular momentum. However, the constraint iVj + Si = even 
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becomes irrelevant in the small U2 limit 1 -. Therefore in this 
limit, the effective low energy Hamiltonian (Eq. [29b can be 
directly mapped to the rotor model (Eq.[T|) with the identifica- 
tion U 2 /2 U, 2N 2 P X /U -> g and /t x Ajf. For the 
biaxial nematic phase to occur, we therefore need a window 
where the conditions AN 2 t 2 x /U2U Q > 1, and Ni x /U < 1 
are simultaneously satisfied. The second condition arises 
here since we need to be away from the superfiuid transi- 
tion point of boson systems. These conditions can be ex- 
pected to be easily satisfied. For example, in a deep lattice 
(V = IOEr where Er is the recoil energy) set by red de- 
tuned light (A = 985 nm) and containing about 10 sodium 
atoms per well, AN 2 t 2 x /U 2 U ~ 14, and Ni x /U ~ 0.4. 

B. Detection of biaxial nematic order 

In this section, we propose a method to detect experimen- 
tally biaxial nematic order in a condensate of spin-one cold 
atoms. Consider the atoms being in a Mott state which has bi- 
axial nematic order parameter. First, as is customary in most 
of the experiments in ultracold atoms 25 , we switch off both the 
lattice potential and the trap, and let the atoms expand freely. 
We then probe the expanding cloud by a right-circular (<r+) 
polarized laser beam. The dielectric tensor of the atoms as 
seen by the laser beam is given by 23 

(ejk) = 8j h + c a (p)d jk - icie jkl (Si) + c 2 (Qjk), (30) 

where the coefficients c a= /o,i,2} depends on the laser fre- 
quency and (p), (S) and (Qjk) are the density, average spin 
and the nematic order parameters of the atomic cloud. If the 
laser frequency is too far detuned from hyperfine splitting fre- 
quency of the atomic levels, C2 vanishes and the nematic order 
is not probed. On the other hand, if the laser frequency is not 
detuned enough from the hyperfine splitting frequency, there 
will be significant absorption which will weaken the intensity 
of the transmitted light. As shown in Ref. |23t there indeed 
exist a window for several spin-one atom species where the 
imaging can be done. 

If the expanding cloud is sufficiently optically thin and ho- 
mogeneous, the polarization of the transmitted beam (taken to 
be propagating along the z axis) is 

pout = e li ^[H / dz (V e (sj/) - 1)] Pin, (31) 

c Jo 

where pj„ = (p x , p y ) is the two-component polarization vec- 
tor of the incoming laser beam, ec xy ) is the reduced dielectric 
tensor in the (xy) plane and A is the thickness of the medium. 
As discussed in Ref. 23, the presence of spin order (S) in 
the atom cloud gives a phase shift to the atoms whereas a ne- 
matic order leads to a left-circular (cr_) polarized component 
in the transmitted beam. So, if we shine on the spin-one sam- 
ple a beam of pure er + light in such a way that the principal 
axes Qi and Q2 of the nematicity ellipsoid are orthogonal to 
the direction of the propagating beam, the intensity of the cr_ 
component of the transmitted beam is given by 

/_ = I a+ ^ £ dz ((Qi) - (Q a » | 2 (32) 



where a + is the amplitude of the incoming beam. Note that 
this method distinguishes between uniaxial and biaxial ne- 
matic ground states. In the uniaxial state {Qi) — (Q2) 
and I- = 0; however, for a biaxial nematic ground state 
(Qi) 7^ {Q2) so J_ 7^ 0. Thus passing the transmitted beam 
through a crossed polarizer, one should be able to measure J_ 
and hence detect the presence of a biaxial nematic state. 



IV. CONCLUSION 

We have studied a disordered 0(2) rotor model with 
quadrupolar interaction and demonstrated that the model ex- 
hibits a biaxial nematic phase in the disordered average sense. 
It is demonstrated that within mean-field analysis, the biax- 
ial nematic phase is stable against small quantum fluctua- 
tions. Such models are shown to be realized in the Mott 
phase of spin-one ultracold bosons in optical lattices with 
spin-dependent disordered potential in the limit of large num- 
ber of bosons per site. We have also suggested an experiment 
which can, using laser imaging of the spin-one atoms, detect 
the biaxial nematic phase. 



Acknowledgments 

This work was supported by NSERC, the Canadian Institute 
for Advanced Research, the Canada Research Chair Program 
(YBK, KS, JSB), and Le Fonds quebecois de la recherche 
sur la nature et les technologies (JSB). YBK thanks Chetan 
Nayak for discussions that sparked his interest in biaxial ne- 
matic phases, KS thanks Duncan O'dell and Ying-Jer Kao for 
helpful discussions. 



APPENDIX A: SPIN-DEPENDENT OPTICAL LATTICE 

We propose here a method to create a spin-dependent opti- 
cal square lattice. Using our approach, trapped bosons with 
S z = { — 1,1} experience the same potential, V—\ = V\, 
while bosons with S z = are subject to a different poten- 
tial Vq. 

Consider atoms with total angular momentum S = 1 in- 
teracting with a configuration of laser beams producing an 
electric field E(r). Building on previous work for 5=1/2 
particles 2 ^, we deduce that atoms with S = 1 experience an 
external potential of the form 

V aP (r) = V{r)5 aP + B(r) • S Q(3 + N l3 {r)G l3aP , (Al) 

with a,0 = {-1, 0, 1}. In Eq.[Al] the scalar potential V(r) 
is proportional to the light intensity, the vector field B(r) is 
proportional to the electromagnetic spin 24 and couples to the 
total atomic angular momentum operator S, and the second- 
rank tensor Ny (r) is proportional to the light nematicity^ 3 - and 
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couples to the quadrupole moment operator Gij : 

V(t) = fe E*(r)-E(r) 
B(r) = -i&!E*(r) x E(r) 

-lE*(r).E(r)5«] 



S Q/3 = (l,«|S|l,/3) 

Gija/3 = 



+ - -S 2 <5, 



(A2) 



The coefficients 60,1,2 are functions of the light frequency and 
the atomic structure. To obtain an effective coupling to the 
light nematicity (i.e. a large enough b 2 value), one needs to 
tune the laser frequency close to the hyperfine splitting of the 
atoms, but far from the fine structure splitting such that b\ <C 

b 2 ,b . 

Then, to generate the optical square lattice (say in the xy 
plane), we use two orthogonal pairs of counter propagating 
monochromatic lasers, and choose these equal intensity light 
fields to be linearly polarized in the z direction. The total 
electric field produced by this configuration is thus given by 



E(t, x, y) = 2 E z e iwt [e l<t >* cos(kx) + e 1 ^ cos(fcy)] , 



where k is the wavevector, and <j> x , <f> y are the initial phases 
for the electric field propagating in the x and y directions re- 
spectively. We choose the difference between these two initial 
phases A<fi = <j> x — <p y to be equal to n/2. Using this elec- 
tric field configuration, we find the electromagnetic spin to be 
zero and the external potential to be 

V a p(x,y) = A\E \ 2 (cos 2 (kx)+ cos 2 (ky)) x 



(b Q --b 2 )5 a p + b 2 (l,a\S 2 \l,l3}). 



(A3) 



(A4) 



Hence, only the diagonal terms of the external potential tensor 
are non-zero and are given by 



V 00 = A(x,y)(b - -62) 
V u = A(x,y)(b + -b 2 ) 
V-1-1 = A(x,y)(b + ^b 2 ), (A5) 



where A(x,y) = 4|_E | 2 (cos 2 (fca;) + cos 2 (fcy)). As a result, 
we obtain a spin-dependent optical square lattice with V\ = 
V-! + V . 
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